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1.  EXECUTIVE  SUMMARY 


This  is  the  final  report  presented  to  the  Survivability  and  Safety  Branch  of  Wright  Laboratory 
(WL/FIVSM).  It  summarizes  efforts  of  the  Satellite  Assessment  Division  of  Phillips  Laboratory 
(PLAVSA)  toward  development  of  a  high  fidelity  predictive  structural  composite  penetration  model. 
This  physics  based  "first  principle"  model  vdll  be  able  to  simulate  the  dynamic  response  and  progressive 
damage  in  structural  composite  laminates  subjected  to  ballistic  impact  of  non-nuclear  treats. 

The  need  for  the  development  of  a  high  fidelity  first  principle  penetration  model  for  composite  structures 
subjected  to  ballistic  impact  is  apparent.  Although  the  smoothed  particle  hydrodynamics  (SPH) 
technique  has  been  successfully  applied  to  simulate  hypervelocity  impact  of  homogeneous  and  isotropic 
solid  materials  in  the  early  1990's,  the  need  for  simulating  behavior  of  non-homogeneous,  anisotropic 
materials  under  intensive  dynamic  loading  remains.  The  new  challenges  posed  by  modeling  the  non¬ 
homogeneity  and  anisotropy  of  fiber-reinforced  composites  are  two-fold.  One  is  how  to  accurately 
describe  anisotropic  elasto-plasticity,  equation  of  state,  and  failure  condition  for  unidirectional 
composites  under  Mgh  strain-rate  loadings.  The  other  is  how  to  feasibly  implement  a  composite  structure 
model  in  a  numerical  code  simulation  while  considering  the  limitation  of  computer  memory  and  time. 

To  address  the  former  problem,  PLAVSA  has  constructed  a  theoretical  foundation  representing  the 
needed  additional  physics  essentially  for  simulating  dynamic  response  of  composite  laminates  under 
severely  elevated  pressure  regimes.  The  remaining  task  is  to  descretize  and  implement  into  the  SPH  the 
governing  equations  where  composite  laminates  can  be  modeled  numerically  feasibly  and  economically. 

Our  results  clearly  demonstrate  success  in  both  accounts:  (a)  numerical  simulation  of  composite 
laminate  damage  zone  correlated  well  qualitatively  with  the  known  experimental  data,  and  (b)  newly 
improved  SPH  proved  to  be  a  robust  and  viable  analytical  tool  for  predicting  response  of  non- 
homogeneous,  anisotropic  composite  materials  to  intensive  dynamic  loadings.  In  conclusion,  we  have 
achieved  the  established  goal  of  this  task  where  a  fundamental  and  economical  analytical  tool  capable 
of  predicting  damage  response  of  fiber-reinforced  composites  to  impact  loadings  is  developed.  We 
submit  however,  there  is  a  great  deal  more  development  needed  before  we  can  assess  damage  across  a 
broad  spectrum  of  composite  structures.  Following  is  a  list  of  most  pertinent  topics: 

(1)  generalization  of  the  current  developed  anisotropic  plasticity  model  and  failure  criterion  for 
fiber-reinforced  composites  to  encompass  high  strain-rate  regime, 

(2)  extension  of  the  proposed  equation  of  state  for  composites  to  include  temperature  and  phase 
changes, 

(3)  development  of  a  tri-axial  ellipsoidal  smoothing  function  for  prismatic  particles,  and 

(4)  extension  of  the  current  SPH  capability  for  interface  and  boundary  condition  applications. 
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2.  INTRODUCTION 


Impact  damage  to  composite  laminates  has  been  studied  extensively  in  the  last  two  decades.  Most 
efforts  have  been  devoted  to  low  velocity  impact.  Due  to  the  constraints  of  technique  and  cost, 
however,  relatively  little  progress  has  been  achieved  in  the  understanding  of  high  velocity  impact 
of  laminated  composites.  Among  the  latter,  Walters  and  Scott  [1]  described  an  experimental 
observation  of  deformation  patterns  within  Kevlar  29®  composite  laminates  subjected  to  a  shaped 
charge  jet.  Brar  [2]  measured  the  force  on  projectile  noses  while  penetrating  composite  targets  as 
well  as  characterized  the  failure  modes  for  S-2  Glass®  composites.  Other  experimental  efforts  focused 
on  the  ballistic  limit  and/or  residual  velocity  of  a  projectile  on  composite  targets  [3-6].  Recently, 
Schonberg  and  Walker  [7]  presented  an  experimental  investigation  of  composite  materials  used  as 
bumper  shields  for  preventing  the  perforation  of  the  pressure  walls  in  multi-wall  structural  systems 
subjected  to  hypervelocity  projectile  impact.  In  the  area  of  analytical  analysis,  Cantwell  and  Morton 
[8]  used  an  energy  dissipation  model  which  accounts  for  elastic  deformation,  delamination  and  shear- 
out  to  predict  the  perforation  threshold  of  a  composite  beam.  This  simple  model  predicted  the 
experimental  trends  successfully  only  for  thin  laminates.  Zhu  et  al.  [9]  developed  a  phenomenonlogical 
model  which  comprised  both  global  and  local  deformations  for  laminated  Kevlar/polyester  plates 
impacted  by  hard-steel  cylindro-conical  projectiles.  Although  an  accord  between  the  prediction  and 
measurement  of  the  ballistic  limits,  residual  velocities  and  displacement  histories  of  the  projectile  was 
satisfactory,  it  deteriorated  successively  when  velocities  and  accelerations  were  examined.  Sun  and  his 
associates  [5,10]  proposed  a  quasi-static  load-displacement  punch  curve  as  the  structural  constitutive 
model  for  the  entire  perforation  process.  The  residual  velocities  were  then  estimated  based  on  the 
balance  of  dynamic  energy.  Vinson  and  Walker  [11]  extended  the  conical  shell  model  for  a  textile 
flexible  cloth  target  [12]  to  predict  the  ballistic  limits  and  residual  velocities  of  the  impactor  for 
fibrous  polymer  matrix  composites.  In  the  area  of  numerical  analysis,  Blanas  [13]  used  the  finite 
element  code  DYNA3D  to  evaluate  ballistic  limits.  He  concluded  that  DYNA3D’s  macromechanical 
model  is  a  marginally  acceptable  method  of  modeling  penetration  phenomena  of  composite  laminates 
only  if  qualitative  results  are  needed  as  a  strucmral  design  criterion.  In  an  attempt  to  model  shock 
wave  propagation  in  fiber-reinforced  composites,  Anderson  and  his  co-workers  [14,15]  modified  the 
Mie-Gruneisen  equation  of  state  (EOS)  for  metals  and  applied  it  to  composites  so  that  it  is  able  to 
account  for  the  portion  of  the  pressure  resulting  from  deviatoric  strains.  The  macromechanics-based 
anisotropic  constitutive  laws  for  composites,  together  with  the  modified  EOS,  were  implemented  into 
the  finite  difference  hydrocode  HEMP  [14]  and  the  finite  element  hydrocode  EPIC  [15].  Since  no  failure 
investigation  was  attempted,  only  the  onset  of  failure  can  be  predicted  in  their  analyses. 

The  phenomena  involved  in  high  velocity  impact  on  composite  laminates  are  complex,  including 
local  contact,  bulging,  global  deformation,  matrix  cracking,  fiber  breakage,  fiber-matrix  interface 

debonding,  delamination,  fragmentation,  etc.  As  pointed  out  by  Sun  and  Potti  [10],  a  detailed 
modeling  of  damage  progression  during  perforation  would  be  a  very  difficult  task.  This  is  why  most 
of  the  above  theoretical  investigations  for  high  velocity  impact  on  composite  materials  was  confined 
to  ballistic  limits  and  residual  velocities  of  a  projectile.  Recently,  Chen  et  al.  [16]  have  applied  the 
SPH  technique  [17]  to  simulate  the  detail  perforation  process  for  B/Al  composites.  Since  the  diameter 
of  boron  fibers  is  relatively  large  as  compared  to  the  lamina  thickness,  they  modeled  the  fibers  and 
matrix  as  two  different  isotropic  materials.  This  approach,  however,  is  prohibitively  computer  memory 
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intensive  for  graphite/epoxy  (Gr/Ep)  composites  because  there  is  a  great  number  of  fibers  over  a  lamina. 

The  new  challenges  posed  by  modeling  the  non-homogeneity  and  anisotropy  of  fiber-reinforced 
composites  are  two-fold.  One  is  how  to  accurately  describe  the  anisotropic  elasto-plasticity,  equation 
of  state  (EOS),  and  failure  conditions  for  a  unidirectional  composite  under  high  strain-rate  loadings.  The 
other  is  how  to  feasibly  implement  a  composite  structure  model  in  a  numerical  code  simulation  while 
considering  the  limitation  of  computer  memory  and  time. 

In  this  report  the  macro-mechanics  approach  for  fiber  composites  was  proposed  for  the  SPH  simulations 
of  impact  penetration  into  Gr/Ep  composites.  A  3-D  elasto-plasticity  constitutive  model,  an  EOS  and 
a  failure  criterion  which  account  for  the  anisotropic  behaviors  of  fiber  composites  were  developed.  Since 
material  often  undergoes  large  rotation  during  impact  process,  stress  and  strain  transformations  between 
the  material  principal  axes  and  the  deformed  configurations  were  also  considered  based  on  the  polar 
stress  rate  approach.  The  above  modules  for  fiber  composites  were  then  incorporated  into  the  SPH. 

3-D  simulations  were  carried  out  with  the  [0/90/45/-45]s  and  [45/-45/0/90]s  Gr/Ep  composite  laminates 
impacted  by  a  steel  cubic  projectile.  The  impact  velocities  exceeding  the  ballistic  limits  were  of  interest. 
Detailed  numerical  results  and  comparison  with  the  experimental  data  will  be  presented  and  discussed. 
Further  development  needs  will  also  be  addressed. 
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3.  THEORETICAL  MODEL  DEVELOPMENTS 

To  describe  the  impact  penetration  behavior  of  homogeneous  solid  media,  three  conservation  laws  of 
mass,  momentum  and  energy  must  be  obeyed.  In  general,  the  conservation  laws  offer  five  scalar 
equations  but  involving  seven  unknowns  which  are  pressure,  density,  internal  energy,  temperature,  and 
three  velocity  components.  In  order  to  obtain  uniqueness  of  the  solutions,  material  properties,  such  as 
constitutive  relationships  and  equations  of  state,  will  have  to  be  imposed.  Since  often  material  would 
undergo  large  rotation  during  the  impact  process,  it  is  also  necessary  to  transform  stresses  and  strains 
back  and  forth  between  the  unrotated  and  deformed  configurations.  In  addition,  a  strength  model  is  ^ 
needed  to  describe  any  possible  fragmentation. 

For  a  fiber  rein-fi)rced  composite,  the  above  physics  are  complicated  and  little  provision  has  been  made 
to  account  for  the  existing  inherent  non-homogeneous  and  anisotropic  nature  so  far.  In  this  section,  a 
generalized,  3-D,  anisotropic  elasto-plasticity  model,  equation  of  state,  and  dynamic  failure  criterion  for 
fiber  composites  will  be  proposed.  The  stress  and  strain  transformations  between  the  material  principal 
coordinates  and  deformed  configurations  will  be  developed. 

3.1  ANTSOTROPTC  ELASTO-PLASTICITY  THEORY 

Several  three-dimensional  anisotropic  plasticity  models  [18-24]  have  been  proposed  for  fiber-reinforced 
composite  materials.  However,  most  of  the  theories  suffer  two  shortcomings.  One  is  that  some 
approaches  utilize  the  plasticity  theory  for  metals  to  describe  the  nonlinear  behavior  of  fiber  composites. 

For  instance,  Grifian  et.  al  [18],  Hansen  et.  al  [19],  and  Xie  and  Adams  [20]  employed  Hill's  orthotropic 
yield  criterion.  It  should  be  noted  that  Hill's  plasticity  theory  assume  that  hydrostatic  stress  does  not 
influence  plastic  deformation  and  that  the  plastic  dilatation  is  incompressible.  These  two  assumptions 
are  particularly  questionable  for  fiber-reinforced  composites.  The  other  shortcoming  is  that  some 
approaches  do  not  describe  material  anisotropy  parameters  completely.  The  assumptions  of  linear 
elasticity  in  the  fiber  direction  [20-22]  and  transverse  isotropy  [20,21,23]  were  commonly  made.  The 
former  is  acceptable  only  for  some  particular  types  of  composites  in  which  the  elastic  fibers  carry  the 
majority  of  load.  The  latter  is  not  suitable  for  composites  with  a  fiber  arrangement  whose  overall 
spacings  are  different  in  the  different  directions.  Recently,  Chen  et.  al  [24]  proposed  a  generalized  yield 
criterion  that  is  quadratic  in  stresses  for  B/Al  and  AS4/PEEK  composites.  This  plasticity  model  not  only 
relaxes  the  previously  imposed  assumptions  [18-23],  but  is  also  general  in  nature  covering  composite 
materials  composed  of  various  fiber  arrays. 

3.1.1  3-D  MTCROMECHANICS  FINITE  ELEMENT  ANALYSIS 

Driven  by  laboratory  testing  limitations  and  cost,  Chen  et  al  [24]  used  the  ANSYS  code  to  conduct  3-D 
micromechanics  finite  element  calculations  to  generate  the  nonlinear  macro  stress-strain  data.  Their 
micromechanics  model  was  based  on  a  representative  volume  element  (RVE)  or  unit  cell.  The  analysis  ^ 
was  performed  under  the  assumption  that  the  periodic  fiber  arrangement  is  in  the  form  of  square  array. 

In  reality,  however,  in  a  composite  lamina  the  fiber  distribution  is  not  regular  within  the  cross-section. 

As  reported  by  Sun  and  Vaidya  [25],  for  the  AS4/PEEK  composite  a  hexagonal  array  gives  a  better 
prediction  of  the  nonlinear  behavior  of  the  composite  than  a  square  array.  In  the  present  3-D 
micromechanics  finite  element  analysis,  thus,  the  hexagonal  array  is  employed  for  the  nonlinear  stress- 
strain  data  acquisition.  The  composite  system  considered  here  is  the  AS4/3501-6  with  66%  fiber  volume 
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ratio  (Fig.  1).  The  fibers  are  assumed  to  be  transversely  isotropic  and  linearly  elastic.  The  matrix 
exhibits  elastic-plastic  characteristics  and  obeys  the  von  Mises  J2  flow  rule.  The  elastic  constants  used 
in  the  numerical  analysis  are  given  as  follows: 


AS4  fiber:  E„  =  234  GPa, 

G23  =  5.5  GPa, 
V23  =  0-25, 


E22  =  E33  =  13.8  GPa, 
Gj2  =  Gi3  =  27.6  GPa, 

Vi2  =  Vi3  =  0.2; 


3501-6  epoxy:  E  =  4.8  GPa, 


V  =  0.34. 


The  nonlinear  uniaxial  stress-strain  curve  for  the  3501-6  epoxy  can  be  found  in  Adams  and  Crane  [26]. 
The  finite  element  model,  load  and  boundary  conditions  applied  to  the  RVE  are  described  in  [24].  The 
numerical  analysis  was  carried  out  using  the  ANSYS  code. 


Figure  2  shows  the  predicted  macro  plastic  strains  in  the  AS4/3501-6  composite  lamina  subjected  to 
hydrostatic  stress.  The  stress-strain  responses  are  the  same  in  both  the  transverse  directions.  It  is  clear 
that  the  sum  of  three  normal  plastic  strain  components  is  not  zero.  This  reveals  that  for  the  AS4/3501-6 
composite  hydrostatic  pressure  not  only  influences  plastic  deformation,  but  the  plastic  dilatation  is 
compressible  as  well.  The  stress-strain  curves  under  uniaxial  simple  tension  loading  in  the  fiber 
direction  are  plotted  in  Fig.  3 .  The  solid  lines  are  the  linearly  elastic  solutions,  and  the  circles  present 
the  elastic-plastic  solutions.  At  1%  of  the  total  longitudinal  strain,  for  example,  the  plastic  strain 
component  is  about  0.35%  of  the  total  strain.  The  contraction  in  the  transverse  directions  at  this  load  is 
about  0.2%.  The  plastic  strain  portion  is  about  7.3%  of  the  total  contraction.  If  the  assumption  of  linear 
elasticity  in  the  fiber  direction  is  posed  in  developing  an  anisotropic  plasticity  model,  then  these  plastic 
deformations  will  vanish.  It  is  recommended,  in  the  interest  of  more  accurate  solution,  that  the  imposed 
assumption  be  removed. 


3.1.2  ANISOTROPIC  YIELD  FUNCTION 


The  generalized,  3-D  anisotropic  yield  function  proposed  by  Chen  et  al.  [24]  is  given  as, 

/(Oy)  =  a, |a,|  + 022^22 ^33*^33’'' 2<3|20],022+ 2^230,2033 
+  2a,30„033+  2a,4023+  2055  03>2agg0f2 

=  k 

where  the  stresses  refer  to  the  principal  material  directions,  and  A:  is  a  state  variable.  The  nine 
plasticity  coefficients  a^j  describe  the  amount  of  anisotropy  in  plasticity.  This  yield  criterion  reduces  to 
the  Hill's  orthotropic  yield  function  when 

^12  =  ^33  -  (0,1 +  ^2+^33)^ 

^13  ”  ^22  ~  (^  1  ]  ^2  ^33^  ^ 

^23  ~  ^11  ~  (^11  ^22^  ^33)^ 

Assuming  the  associated  flow  rule,  the  incremental  plastic  strains  Jcy  can  be  written  as 
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Fig.  3 


(3) 


dd  =  dk 


So.., 


in  which  the  superscript  p  denotes  plasticity,  and  dX  is  a  plastic  factor.  Let  the  effective  stress  a  be 
defined  as 


a  = 


N  2 


•/ 


By  comparing  Eqs.  (1)  and  (4),  one  finds  that 


7  2—2 

k  =  —a 


(4) 


(5) 


The  incremental  effective  plastic  strain  can  be  derived  through  the  concept  of  plastic  work: 


dW’’  -  Oydey  =  add’ 


(6) 


From  Eqs.  (1),  (3)  and  (6),  thus,  the  incremental  effective  plastic  strain  can  be  expressed  as 

de^  =  —adk 


(7) 


Substituting  Eq.  (1)  into  (7)  with  the  inversion  of  Eq.  (3)  for  Oy  and  regrouping  it  yields  de^  in  an 
explicit  form 


(de^^  =  -^^Aj,(def,f+A22(de^2)^+Aj^(de^3f+  2A^^deF^^de^^+  lA^^dd^^de^^^ 


(8) 


where 


Oj,  0,2  0,3 


^13  “^23  ^33 


(9) 


and 


2  2  2 
.(4,,— O22O33- O23,  .<422“  0,,033  -  0,3  ,  -^33””  422**^  12 

-^12“  ^13‘^23’  ‘^12^33’  ^13 “  *^12^23  '  ^13*^22  ’  ^23“  *  ®11^3 


(10) 


^  The  determinant  in  Eq.  (9)  is  zero  for  a  material  which  obeys  the  von  Mises  flow  rule  or  the  Hill 

criterion  because  of  the  condition  of  the  incompressibility  of  plastic  dilatation. 

The  o-  d’  relationship  for  a  particular  loading  can  then  be  obtained  from  Eqs.  (1),  (3),  (4)  and  (8).  For 
example,  consider  a  uniaxial  loading,  a,,  5^0  only  (no  summation).  By  using  Eqs.  (1)  and  (4),  the 
effective  stress  becomes 
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a  = 


3(3.. 


N  2 


(11) 


The  incremental  plastic  strain  components,  Eq.  (3),  can  be  written  as 

a. 


(k!'-  = 


del,  del=^del  del  =  0 


Substitution  of  Eq.  (12)  into  (8)  yields  the  incremental  effective  plastic  strain, 


(12) 


(13) 


The  incremental  eflFective  plastic  strain  can  also  be  derived  by  substituting  (1 1)  into  (7)  and  utilizing  (3). 


For  a  shear  loading,  0^*0  only,  the  effective  stress  becomes 


(14) 


where  rr  =  44,  55  or  66  depending  upon  the  shear  stress  component  (see  Eq.  1).  The  incremental 
effective  plastic  strain  is  then  given  as 


deP  = 


2deP 


(15) 


In  principle,  a  master  o-e^  relationship  for  a  unidirectional  composite  can  be  defined  with  any  of  the 
three  normal  stress-strain  curves  and  three  shear  stress-strain  data,  using  either  Eqs.  (11)  and  (13)  or 
Eqs.  (14)  and  (15).  In  this  study,  the  macro  stress-strain  response  in  the  2-^rection  together  with  settii^ 
^22=  1.0,  without  loss  of  generality,  was  considered  to  define  the  master  o-e^  curve.  Once  the  master  o-e^ 
curve  for  the  composite  is  obtained,  the  values  of  a^i,  a^,  and  can  be  chosen  by  the  trial  and 

error  opfimization  such  that  all  five  stress-strain  curves  are  brought  into  coincidence  with  the 
master  o-  curve.  The  other  constants  associated  with  the  interaction  between  any  two  normal  stresses, 
ai2,  ^13  and  ^23,  can  be  approximated  using  Eq.  (12).  The  detail  of  the  calibration  procedure  can  be  found 
in  Ref  9. 

Using  the  predicted  macro  stress-strain  data  from  the  micromechanics  finite  element  analysis,  the  nine 
plasticity  parameters  in  the  anisotropic  yield  function  (Eq.  (1))  for  the  AS4/3501-6  composites  are 
calibrated  as  follows: 


=  0.00058, 

^22  ^33 

0.975, 

ai2  =  ai3  =  -0.0025 

=  0.58, 

^55  ~  ^66  ”  1-45. 

> 

The  corresponding  master  effective  stress-effective  plastic  strain  curve  is  shown  in  Fig.  4. 

3.1.3  VERTFTCATTON  OF  THE  PLASTICITY  MODEL 

The  verification  of  the  anisotropic  plasticity  model  established  for  the  AS4/3501-6  composite  was 
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Fig.  4  A  master  effective  stress-effective  plastic  strain  curve 


conducted  with  the  three  cases:  (1)  uniaxial  tension  in  the  2-direction,  (2)  uniform  dilatation,  and  (3) 
a  tri-axial  stress  case  with  the  loading  ratio  10:l:-2.  The  effective  moduli  used  in  the  verifications  are 
predicted  from  the  3-D  micromechanics  finite  element  model  and  are  given  as  follows; 

E,i  =  1 56  GPa,  E22  =  E33  =  9.28  GPa, 

G23  =  3 .78  GPa,  Gi2  =  G13  =  7.36  GPa, 

V23  —  0.389,  Vj2  =  V 13  =  0.244. 

Figures  5  shows  the  comparison  of  the  total  stress  and  strain  results  predicted  by  the  3-D 
micromechanics  finite  element  model  and  the  anisotropic  plasticity  mode  for  the  uniaxial  tension  in  the 
2-direction.  The  agreement  between  the  two  results  is  excellent.  At  022=  0.209  GPa,  the  relative  errors 
are  1 .6%  for  622  and  0.9%  for  .  The  uniform  dilatation  case  is  presented  in  Fig.  6.  As  seen  in  the 
figure,  the  solutions  predicted  by  the  plasticity  model  match  very  well  with  those  predicted  by  the  3-D 
micromechanics  finite  element  model.  It  is  interesting  to  note  that  the  plasticity  in  this  case  is 
insignificant.  Figures  7  shows  the  comparison  of  the  total  stress  and  strain  results  for  the  tri-axial  stress 
case.  Again,  there  is  an  excellent  agreement  between  both  predictions.  At  the  maximum  load  calculated 
here,  the  relative  errors  are  2.4%  for  e,, ,  2. 1%  for  €22 ,  and  1 . 1%  for  €33 . 

3.1.4  ET.ASTO-PLASTTC  TANGENT  STIFFNESS  MATRIX 


With  the  confidence  of  the  anisotropic  plasticity  model,  the  elastic-plastic  constitutive  relationship  will 
be  established.  The  plastic  factor  dX  can  be  determined  by  rewriting  Eq.  (7) 

dX  =  1-^  (16) 

in  which  =  do/dd’  is  the  slope  of  the  master  a-e^  curve. 


If  the  incremental  strains  are  small,  one  is  able  to  linearly  decompose  them  into  the  elastic  part  d€^ 

and  the  plastic  part  de^  as 

The  elastic  strain  increment  is  defined  by  the  linearly  elastic  constitutive  relationship 

{d^}  = 

where  [5'']  is  the  elastic  compliance  matrix,  and  {do}  is  the  incremental  stress  vector.  From  Eqs.  (3)  and 
(16),  one  has  the  relations  between  the  incremental  plastic  strains  and  the  incremental  stresses  as 

{deP}  ^  [SP]{do} 


where  [5 'I  is  the  plastic  compliance  matrix.  The  entities  of  [5^]  are  obtained  by  substituting  Eqs.  (1), 
(4)  and  (16)  for /,  do,  and  dX,  respectively,  into  Eq.  (3).  They  are 


in  which 


vP 


pC,C 


i,j=  1,2,. ..,6 


do.  ’  ^  4  3^ 


(20) 


(21) 


12 
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Stress-strain  results  under  simple  tension  in  the  transverse  direction  predicted  by  the  micromechanics 
finite  element  model  and  plasticity  model 
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Fig.  6  Stress-strain  results  under  uniform  dilatation  predicted  by  the  micromechanics 
finite  element  model  and  plasticity  model 
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Stress-strain  results  under  tri-axial  loading  predicted  by  the  micromechanics 
finite  element  model  and  plasticity  model 


where  the  subscript  /  for  o:  1=11,  2=22,  3=33,  4=23,  5=31,  and  6=12.  The  relations  of  the  complete 
elastic-plastic  stress  increments  c?o.  and  the  total  strain  increments  c?e..  are  then  derived: 

{do}  =  {D^P]{de}  ^22) 

where  the  elastic-plastic  tangent  stiffness  matrix  is 

[D^P]  =  ([5"]  +  [SP])-^ 

Equation  (22)  can  uniquely  determine  the  elasto-plastic  stress  changes  {da}  which  arise  during  any 
iteration  in  which  known,  finite  changes  of  strain  {de}  are  imposed. 

To  calculate  the  elasto-plastic  stresses,  the  scaling  factor  approach  [27]  is  employed.  At  each  time  step, 
atrial  stress  state  is  first  computed  with  the  assumption  of  elastic  deformation.  If  the  loading  surface  is 
violated  by  the  elastic  trial  stresses,  then  the  stresses  will  be  bring  back  to  the  yield  surface  from  the 
overshooting  stress  state.  In  this  study,  the  portion  of  the  strain  increment  corresponding  to  the  stress 
increment  beyond  the  yield  value  is  divided  into  5  sub-intervals.  The  final  stresses  are  determined  by 
summing  the  stresses  on  the  yield  surface  and  the  elasto-plastic  stress  increments  sequentially 
calculated  five  times,  using  Eq.  (22)  with  the  updated  plastic  compliance  matrix  [D^p]  each  time. 

3.2  FOTIATTON  OF  STATE 

For  metal  materials,  plastic  flow  is  often  assumed  to  not  be  influenced  by  hydrostatic  pressure.  Thus, 
stresses  are  split  into  two  parts,  i.e.,  hydrostatic  pressure  and  deviatoric  stresses.  The  hydrostatic 
pressure/?  is  then  calculated  using  an  equation  of  state  (EOS)  having  the  functional  form/?  — /?(p ,  e). 

For  a  fiber-reinforced  composite,  both  the  classical  assumptions  used  in  the  Hill's  orthotropic  yield 
function  are  questioned  as  shown  in  Fig.  2.  Thus,  the  application  of  the  conventional  approach  in  which 
the  total  stresses  are  partitioned  into  hydrostatic  pressure  and  deviatoric  stresses  to  fiber-reinforced 
composites  remains  a  challenge. 

Recently,  Anderson  et.  al  [15]  proposed  a  modified  version  of  the  Mie-Gruneisen  EOS  for  fiber 
composites.  The  pressure  used  in  their  EOS  is  still  defined  as  the  negative  of  the  average  of  the  three 
normal  stress  components.  Therefore,  the  modified  EOS  must  account  for  the  portion  of  the  pressure 
resulting  from  the  deviatoric  strains,  which  is  evaluated  based  on  the  linear  elasticity  assumption.  The 
credibility  of  this  EOS  was  not  validated  yet,  especially  the  linear  elasticity  correction  for  the  deviatoric 
strains. 

It  is  interesting  to  notice  the  results  shown  in  Fig.  6  for  the  AS4/3501-6  composite  under  uniform 
dilatation.  Although  some  plasticity  exists  in  the  composite,  it  is  much  less  pronounced  as  compared 
with  that  under  hydrostatic  stress  (Fig.  2).  Therefore,  the  assumption  that  uniform  dilatation  does  not 
influence  plastic  deformation  a  better  approximation  for  the  AS4/3501-6  composites.  With  this 
assumption,  it  may  be  convenient  to  define  the  "pressure"  as  the  negative  of  the  average  of  the  three 
normal  stress  components  only  corresponding  to  the  volumetric  strain.  Then  this  pressure  can  be 
computed  using  the  Mie-Gruneisen  EOS  without  any  modifications.  Once  the  pressure  change 
between  two  consecutive  time  steps  is  determined,  the  three  normal  stress  components  corresponding 
to  the  change  of  uniform  dilatation  can  be  retrieved  using  the  linearly  elastic  relationship.  The  total  trial 
stress  increments  are  then  evaluated  by  summing  the  normal  stress  changes  and  the  elastic  "deviatoric" 
stress  increments  which  are  associated  only  with  the  deviatoric  strains.  Then  the  scaling  factor  algorithm 
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mentioned  previously  is  employed  to  compute  the  elasto-plastic  stresses.  Note  that  the  "pressure' 
introduced  here  is  different  from  the  hydrostatic  pressure  defined  by  Anderson  et  al  [15], 


The  Me-Gruneisen  EOS  is  a  theoretically  sound  model  for  crystalline  metals.  It  is  suitable  for  high  (but 
not  extremely  high)  velocity  impact  events  [28],  The  Mie-Gruneisen  EOS  is  given  as: 

p(p,e)  =  (1  - +  Tpe  (24) 

where  T  is  the  Gruneisen  constant  of  the  material,  r|  =  p/po-1  represents  the  compressive  behavior,  and 
/7h(p)  is  the  so-called  Hugoniot  curve  which  represents  the  locus  of  all  points  that  may  be  reached  by 
a  shock  transition  from  an  initial  unshocked  state  (Po,Po)-  Several  forms  of  the  Hugoniot  have  been 
proposed.  The  one  commonly  used  in  SPH  is  written  as 

^  LjTI  +  +  Sh'  h>0  (25) 

fljTl  T|  <  0 

The  coefficients  and  in  Eq.  (25)  are  obtained  from  experimental  shock  compression  data. 

Alternatively,  they  are  analytically  determined  through  a  Taylor's  series  expansion  of  the  Hugoniot 
curve  Pij  =  c^T[{\  +  r\)l[\-'  ti(S-  1)^]  as  follows  [17]: 

=  a[{\+2{l-\)]  (26) 

=  «,[2(M) +3(M)2] 


in  which  ^  is  the  constant  in  the  approximated  linear  relationship  between  the  shock  velocity  (mJ  and 
particle  velocity  (w^). 


=  r  +  (27) 

Since  the  sound  speeds  are  different  in  the  three  directions  for  fiber  composites,  an  equivalent  (bulk) 
sound  speed  c  used  in  Eqs.  (26)  and  (27)  will  be  proposed.  For  a  fiber  composite,  the  linear  bulk 
modulus  K  evaluated  based  on  the  uniform  dilatation  condition  is  given  as 


K  = 


1  , 

l-'’23'’32  +  1-V,3v31  ^  l-V.jVj,  ^  ^ 

\ 

'’21+ ''31 '’23  ^  '’31+ '’21 '’32  ^  '’ 12+ '’ 12 '’3I 

9A 

£1-  £, ,  £„  E. , 

22  33  1 1  33  1 1  22 

^  £22^33  ^22^33  ^11^11  ; 

(28) 


Equation  (28)  can  be  simplified  for  isotropic  materials  as, 

^  E 
3(l-2v) 


(29) 


The  equivalent  (bulk)  Poisson's  ratio  for  the  composite  is  assumed  to  be 


V 


v.(l-F,) 


(30) 


Thus,  the  equivalent  sound  speed  c  for  the  composite  can  be  determined  using  the  relationship  from 
isotropic  materials 


c  = 


£(i-v)_  _ 

(1+  v)(l-2v)p 


(31) 


where  p  is  the  density  of  the  composite,  and  E  is  the  equivalent  Young's  modulus  of  the  composite  and 
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is  evaluated  using  Eq.  (29)  with  the  bulk  modulus  and  equivalent  Poisson's  ratio. 

The  constants,  F  and  and  the  bulk  sound  speed  for  the  AS4/3506-1  composite  used  in  the  EOS 
calculation  are  set  as  follows:  F  =  0.87,  ^  =  1.57,  and  c=  4.69 km/s. 

It  is  important  to  note  that  the  proposed  Mie-Gruneisen  EOS  for  anisotropic  materials  as  well  as  the 
determination  of  the  constants  a,,  a2,  and  flj  using  Eq.  (26)  should  be  validated  experimentally.  In  this 
study,  the  "pressure"  was  computed  using  the  EOS  for  updating  the  shock  speed.  It  was  not  used  in  the 
stress  calculation. 


3.3.  DYNAMIC  FAILURE 


The  theoretical  investigation  for  detail  dynamic  damage  in  composite  structures  due  to  high  velocity 
impact  is  very  limited.  One  of  the  main  reasons  is  lack  of  a  high  fidelity  rate-dependent  failure  criterion. 
Recently,  Randles  and  Memes  [29]  proposed  a  continuum  damage  model  for  thick  quasi-isotropic 
laminates  based  on  the  assumption  of  homogenous,  transversely  isotropic  behavior  for  both  stiffhess  and 
damage.  They  also  demonstrated  a  good  prediction  of  delamination  type  damage  and  spallation  under 
uniaxial  strain  condition  [30].  It  should  be  noted  that  the  approximation  of  transverse  isotropy  is  suitable 
only  for  stiflfiiess.  Once  local  damage  occurs,  the  inital  quasi-isotropic  laminate  will  lose  its  transversely 
isotropy.  Therefore,  the  capability  of  their  continuum  damage  model  for  detailed,  3-D  damage 
predictions  should  further  be  verified. 


To  predict  detailed  damage  in  a  fiber  composite,  using  a  layer-wise  damage  model  in  the  SPH  could  be 
a  more  practical  approach.  Since  it  is  still  not  available  so  far,  a  modified  maximum  stress  failure 
criterion  is  employed  for  the  dynamic  failure  analysis  in  this  study.  Fracture  of  a  particle  is  said  to  have 
occurred  if  any  tensile  stress  or  shear  stress  in  principal  material  directions  is  greater  than  or  equal  to 
the  corresponding  strength. 


Ojj 


12l 


'131 


231 


^  S 

>  s 

>  s. 


12 


13 


23 


(32) 


where  X,  and  Y,  are  the  tensile  strengths  in  the  fiber  and  transverse  directions,  respectively,  and  are 
the  shear  strengths.  When  any  one  of  the  inequalities  in  Eq.  (32)  is  satisfied,  the  particle  is  assumed  to 
instantaneous  fail  in  the  respective  mode.  Once  it  happens,  the  material  can  no  longer  support  the 
corresponding  tensile  and  shear  loads.  For  instance,  if  then  o,,=  0,2=  0,3=  0.  If  I ‘^23  1  >^23  >  then 

023=  0,  022=0  if  O22>0  ,  and  033=0  if  o^^>0.  No  compressive  failure  is  assumed. 

In  reality,  dynamic  fi-acture  strengths  depend  on  strain  rate,  temperature,  and/or  other  factors.  Therefore, 
a  range  of  rate-dependent  stress-strain  curves  would  be  necessary  to  apply  for  simulating  a  high  rate 
interaction  event.  However,  it  is  difficult,  technically,  to  obtain  a  full  range  of  temperature-  and  rate- 
dependent  material  properties,  especially  for  the  strain  rate  beyond  10^  cm/cm-s.  For  simplicity,  only 
one  value  was  assigned  for  each  strength  in  this  study.  Because  of  the  ductile  behavior  in  the  epoxy 
resin,  the  dynamic  transverse  strength  and  shear  strengths  of  the  composite  were  assumed  to  be  five 
times  the  static  strengths.  On  the  other  hand,  the  dynamic  longitudinal  strength  remained  its  static  value 
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since  the  plasticity  in  the  fiber  direction  is  much  insignificant.  Thus,  the  values  of  the  dynamic  strengths 
for  the  AS4/3501-6  composite  used  in  this  study  are:  X,=  \  .447  GPa,  7,  =  0.26  GPa,  and  S,2  =  Sj3  =  S23 
=  0.465  GPa. 


3.4  STRESS  AND  STRAIN  TRANSFORMATIONS 


In  the  analysis  of  composite  laminate  penetration,  factors  which  contribute  to  misalignment  of  the  global 
geometry  coordinate  system  (GGCS)  and  material  coordinate  system  (MCS)  include  initial  layups  and 
large  rotation  which  material  experiences  during  the  penetration  process.  It  is  convenient  to  evaluate 
stresses  in  the  MCS  and  then  transform  them  back  to  the  GGCS  for  displacement  and  subsequent  strain 
calculations.  Therefore,  it  is  necessary  to  trace  the  principal  material  axes  and  transform  stresses  and 
strains  back  and  forth  between  the  GGCS  and  MCS,  during  the  entire  solution  course.  In  doing  that,  the 
Green-Naghdi  stress  rate  approach  was  adopted  because  it  is  supperior  to  the  Jaumann  stress  rate 
approach.  Our  transformation  algorithm  for  fiber  composites  is  extended  fi'om  the  scheme  of  Flanagan 
and  Taylor  [3 1]  which  was  developed  for  isotropic  materials. 


Let  the  GGCS  and  MCS  be  two  sets  of  rectangular  Cartesian  coordinates.  The  rotation  between  these 
two  systems  can  be  specified  by  the  three  Euler  angles  (Fig.  8).  The  individual  transformation  matrices 
for  the  three  consecutive  rotations  are  anti-symmetric  and  can  easily  be  derived.  The  transformation 
matrix  for  the  initial  material  layups,  which  is  from  the  original  (fixed  reference)  GGCS  (Xo,yo,Zo) 
to  the  MCS  (x3,y3,Z3),  is  the  chain-product  of  the  three  individual  transformation  matrices  and  is  given 
as 


CjCj  Cji'j  ■SiCjCj  5,53 


5j  C3  +  C,  ^  l‘*2 


.^2^3 


(33) 


where  c.  =  cos0. ,  s.  =  sin0. ,  and  0,  is  the  /-th  Euler  angle. 

Let  [F]  be  the  deformation  gradient.  Applying  the  polar  decomposition  theorem  to  [F]  yields 

[F]  =  [F][F]^,  =  [F]„,[t/]  (34) 

where  [V]  and  {U\  are  the  left-  and  right-symmetric  stretch  tensors,  respectively,  and  is  an 
orthogonal,  material  rotation  tensor.  The  second  part  of  Eq.  (34),  [F]o,.[C7],  may  be  viewed  as 
deformation  being  first  purely  stretched  from  the  reference  configuration  to  an  unrotated  configuration 
and  then  purely  rotated  from  the  unrotated  configuration  to  the  current  configuration.  The  first  part, 
[V][i?]„^,  may  be  viewed  as  deformation  being  rotated  from  the  reference  configuration  to  an  rotated 
configuration  and  then  stretched  from  the  rotated  configuration  to  the  current  configuration.  Thus,  [f/j 
is  the  stretch  referred  to  the  reference  GGCS  and  [V]  the  rotated  frame. 

The  velocity  gradient  denoted  by  [L]  may  be  expressed  in  terms  of  [F]  as 

[Z]  =  [F][F]->  (35) 

where  the  superposed  dot  indicates  the  time  derivative  and  the  superscript  "-1"  denotes  the  inverse  of 
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the  tensor.  Conventionally,  the  velocity  gradient  is  simply  partitioned  into  the  strain-rate  (symmetric) 
tensor  [D]  and  spin  (anti-symmetric)  tensor  [W], 

[Z]  =  [D]  +  [W\  (36) 

It  is  noted  from  the  definition  of  [Z]  that  [D]  and  [W]  depend  on  the  instantaneous  rate  of  deformation 
of  the  current  configuration.  On  the  other  hand,  [Z]  can  also  be  expressed  in  terms  of  [f/]  and  [i?]  which 
relate  the  current  configuration  to  the  reference  GGCS  by  using  the  right  decomposition  from  (34)  in 
(35), 


[L]  =  +  [Q] 


(37) 


where 


[LI  =  [U][U]'  and  [Q]  =  [/?]„,  [i?];,  (38) 

Substitution  of  the  identity  matrix,  [/],  for  [i?]  into  Eqs.  (37)  and  (38)  indicates  that  [ZJ^  is  a  velocity 
gradient  referred  to  the  unrotated  (not  reference)  configuration.  Thus,  similar  to  Eq.  (36),  [ZJ^  may  be 
partitioned  into  the  strain  rate  tensor  [D\„  and  spin  tensor  [W]o  in  the  unrotated  frame 

[LI  =  [Dl  +  [W\o  (39) 

Also,  the  following  relationships  can  be  derived 


[D]-- 


and  m-[W]-  [R\AW]„mo 


(40) 


It  can  be  seen  from  Eq.  (37)  by  letteing  [C/]=[/]  in  Eq.  (38)  that  [Q]  represents  the  rate  of  rigid-body 
rotation  at  a  material  point.  It  is  equally  simple  to  show  that  [W]  represents  the  rate  of  rotation  of  the 
principal  axes  of  the  strain  rate  [D],  Both  rotation  tensors  are  anti-symmetric.  However,  [Q]  =  [W]  if  and 
only  if  [ZJo  is  symmetric,  which  in  general  fails  to  be  so  although  [17]  is  always  symmetric. 


With  [D]  and  [W],  the  objective  Jaumann  stress  rate  [6]  is  written  as. 


and 


[61=  [dl-[fV\[o]+[ol[m 


[6],=  T  ([oL[Z)]) 


(41) 

(42) 


where  [o]^  is  the  Cauchy  stress  in  the  current  configuration  and  in  general,  is  an  anisotropic  function 
of  a  symmetric  2nd-order  tensor. 

As  proved  by  Johnson  and  Bammann  [32],  the  Cauchy  stress  rate  referred  to  the  unrotated  configuration 
is  also  material  frame  objective.  The  relationships  among  the  unrotated  and  rotated  Cauchy  stress  rates,  [d]^^ 
and  [6]^,  and  the  Green-Naghdi  stress  rate  [b]^  are  given  as 

[ol=  [b],-[Q][o]^+[o]JQ]  =  [7?]^[d]Ji?t 


The  two  Cauchy  stresses  are  related  by 


[o]„=  mUoURl 


(44) 


Note  that  the  Green-Naghdi  stress  rate  and  the  Jaumann  rate  are  very  similar  in  form.  Although  they  are 
all  objective,  the  advantage  of  the  unrotated  Cauchy  stress  rate  and  the  Green-Naghdi  rate  over  the 
Jaumann  rate  is  that  the  first  two  are  truly  a  measure  of  the  rate  of  change  of  stress,  while  the  last  has 
no  conjugate  measure  of  finite  strain  associated  with  it.  In  addition,  the  Jaumann  rate  will  lead  to  a 
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nonsymmetric  stifl&iess  tensor  between  the  rates  of  stress  and  deformation  unless  the  material  is 
incompressible  [32], 


The  rate  type  constitutive  equations  for  the  unrotated  Cauchy  stress  rate  and  Green-Naghdi  rate  can  be 
generally  expressed  as 

K=  T^([o],,  [£>],)  and  [6]^  =  T^([o]^,[Z)]) 


The  second  constitutive  law  in  Eq.  (45)  is  identical  to  Eq.  (42);  however,  the  rotated  Cauchy  stress  rate 
calculated  from  Eqs.  (41)  and  (43)  is,  in  general,  different.  Equation  (45)  can  be  used  to  evaluate  stress 
rate,  but  it  is  convenient  to  evaluate  the  stresses  and  then  examine  the  failure  modes  in  the  material 
frame  for  fiber  composites.  This  can  be  achieved  by  rotating  the  strain-rate  [D]  and  the  previous  Cauchy 
stress  tensors  back  to  the  MCS.  The  transformation  matrix  [R]  from  the  current  frame  to  the  MCS  is 

simply  the  product  of  the  two  matrices 


[i?]  =  [RURl 


(46) 


in  which  is  given  in  Eq.  (33).  Of  consequence,  the  following  transformations  exist 

[o  ]„  =  [Rf[o  ],[/?]  and  [a  =  [R]  [o  ]„ 

[D]^  =  and  [£>]  = 


The  subscript  m  in  the  above  equations  represents  the  quantities  in  the  material  frame.  With  [!)]„ 
obtained  from  Eq.  (48),  the  increments  of  volumetric  and  deviatoric  strains  in  the  material  coordinates 
ran  be  decomposed.  Therefore,  one  can  calculate  the  stress  rate  using  the  EOS  and  the  anisotropic 
elasto-plasticity  theory  developed  in  Section  2.1,  and  then  update  the  current  stress  [o],„ .  As  pointed 
out  by  Flanagan  and  Taylor  [31],  one  of  the  most  challenging  tasks  in  large  deformation  analysis  is  to 
determine  the  rotation  tensor  Multiplying  the  second  equation  in  Eq.  (38)  by  yields, 

[4.= 


Direct  integration  of  Eq.  (49)  gives  the  total  rotation  tensor  [7?]^.  However,  the  orthogonality  of  [R]^ 
degenerates  rapidly  no  matter  how  fine  a  time  increment  used  in  the  integration  [31]. 

In  the  SPH  numerical  analysis,  field  variables  are  evolved  incrementally.  The  geometry  configuration 
is  updated  every  time  step,  while  the  reference  CjGCS  remains  the  same.  Thus,  the  above  transformation 
algorithm  can  be  considered  as  for  a  single  time  step.  To  trace  the  principal  material  axes,  the  current, 
total  rotation  [R]^  at  each  SPH  particle  should  be  estimated  during  the  entire  solution  course. 


An  alternative  approximation  algorithm  for  the  rate  integration  of  [i?]^,:,  which  was  proposed  by  Hughes 
and  Winget  [33],  was  employed.  The  total  rotation  [/?]^  at  time  t+At  is  updated  via 

IK  Jo.  -  [Qjn.  (50) 

where  =  [0J  {x,} .  With  the  assumption  of  a  constant  rate  of  rotation  over  a  time  increment,  the 
approximation  of  Eq.  (50)  is  given  as 

([/]  -  ^Ar[Q])[/?,,^,]„,  =  ([/]  +  iA/[Q])[7?^,  (51) 

The  initial  conditions  for  [i?]„^  is  the  unit  identity  matrix. 
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Before  the  calculation  of  Eq.  (51),  [Q]  must  be  determined.  Substituting  Eq.  (34)  into  (35), 
postmultiplying  by  [V],  and  regrouping  results  in 

[F]=  [Z][F]  -  [F][Q]  (52) 


Since  the  tensors  [Q]  and  [1^  have  only  three  independent  components,  for  convenience,  two  dual 
vectors  based  on  the  right-handed  rates  of  rotation  can  be  defined 

Qy  =  (53) 

W.J  =  (54) 

Substituting  Eqs.  (36),  (53)  and  (54)  into  (52)  results  in 


{0)}  =  {w}  +  mtr{[V\)  -  [F])-'{2} 

in  which 

z.  =  e  .  D  V  , 

/  ikj  jm  mk 


(55) 

(56) 


Making  use  of  Eqs.  (53),  (55)  and  (56),  the  rotation  rate  tensor  [Q]  can  be  evaluated  once  the  velocity 
gradient  [L]  and  the  left  stretch  tensor  [F]  are  determined.  Note  that  the  matrix  [L]  is  computed  using 
the  velocity  components  obtained  fi-om  the  momentum  equations.  The  matrix  [F]  is  time  integrated  using 
Eq.  (52)  with  an  initial  condition  corresponding  to  the  initial  state  of  stress.  For  an  undeformed  state, 
normally  the  stretch  is  set  to  the  identity  tensor  [7]. 


3.5  SMOOTHED  PARTICLE  HYDRODYNAMICS 

SPH  is  a  pure  Lagrangian  particle  method,  which  was  introduced  by  Lucy  [34]  and  (jingold  and 
Monaghon  [35]  in  1977.  It  has  been  successfully  applied  to  astrophysics  and  shock  physics  since  then 
as  well  as  to  hypervelocity  impact  of  solid  materials  in  the  early  1990's  [17].  The  power  of  SPH  is  that 
no  underlying  grids  are  required.  It  not  only  avoids  the  numerical  difficulties  of  mesh  entanglement  and 
distortion  during  the  loading  process  as  often  occur  in  the  finite  element  and  finite  difference  analyses 
but  also  lends  itself  to  the  treatment  of  initiation  and  growth  of  cracks.  Thus,  in  addition  to  the  ballistic 
limits  and  residual  velocities,  the  detailed  damage  progression,  such  as  perforation,  matrix  cracking, 
fiber  breakage,  delamination,  fragmentation,  etc.,  can  naturally  be  simulated  with  this  unique  technique. 

The  foundation  of  SPH  is  the  interpolation  theory.  Through  the  use  of  kernel  estimates,  the  partial 
differential  equations  of  the  three  conservative  laws  of  continuum  mechanics  can  be  transformed  into 
integral  forms.  Numerically,  these  integral  equations  are  approximated  in  terms  of  the  field  variables 
at  a  set  of  disordered  points  (particles).  These  field  variables  are  then  directly  evolved  from  the 
interactions  among  particles  using  an  explicit  time  integration  scheme.  Therefore,  no  matrix  equation 
solving  is  required.  The  details  of  the  SPH  formulation  for  solid  materials  can  be  found  in  Libersky  et 
al.  [17]. 


In  the  present  study,  the  capability  of  the  SPH  for  isotropic  materials  is  extended  to  include  fiber- 
reinforced  composites.  The  anisotropic,  elasto-plastic  constitutive  law,  equation  of  state,  rate-dependent 
failure  criterion  for  the  AS4/3501-6  composites,  as  well  as  the  stress  and  strain  transformations  between 
material  and  geometric  coordinate  systems  under  large  rotation  described  in  the  previous  sections  were 
incorporated  into  the  hydro  code  MAGI. 
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4.  NUMERICAL  RESULTS  AND  DISCUSSION 


The  main  puqjose  of  this  program  is  study  the  feasibility  of  the  SPH  for  3-D  composite  laminate 
penetration.  Before  doing  that,  the  credibility  of  the  modified  MAGI  code,  named  CMAGI,  was  first 
verified.  The  resolution  dependency  of  the  number  of  particles  on  the  SPH  solution  was  also 
investigated.  An  aluminum  plate  impacted  by  a  steel  rod  was  considered  for  the  former,  and  the 
AS4/3501-6  composite  laminates  impacted  by  a  steel  rod  were  studied  for  the  latter.  These  calculations 
were  performed  under  the  plane  strain  assumption.  In  the  3-D  simulations,  normal  impact  of  a  steel  cube 
striking  at  the  center  of  the  8-layer  AS4/3501-6  composite  laminates  were  conducted.  Three  impact 
velocities  (V)  were  investigated:  0.2  km/s,  0.6  km/s  and  2.0  km/s.  For  the  sake  of  clarity,  the  results  of 
damage  in  the  targets  presented  in  the  following  are  confined  to  specific  critical  regions. 

4.1  VATTPATTON  OF  THE  MODIFIED  MAGI  CODE 

The  credibility  of  the  CMAGI  code  for  the  composite  materials  was  verified  with  the  case  of  normal 
impact  of  a  square  steel  rod  striking  at  the  center  of  an  aluminum  plate.  Plane  strain  condition  was 
assumed.  The  width  and  thickness  of  the  cross-section  of  the  target  are  1.93  cm  and  0.0965  cm, 
respectively.  The  size  of  the  square  projectile  is  0.386  cm.  The  impact  speed  was  set  at  2.0  km/s.  Due 
to  symmetry,  only  half  of  the  domain  was  analyzed.  Two  calculations  were  carried  out  using  the  MAGI 
code  (isotropic  material  version)  and  the  CMAGI  code,  respectively.  The  aluminum  target  was  modeled 
as  a  single  material  for  the  former  while  it  was  treated  as  an  8-layer  "composite"  for  the  latter.  For  the 
comparison  purpose,  no  failure  was  assumed  for  the  aluminum  target  since  the  failure  model  proposed 
in  the  CMAGI  is  different  from  that  used  in  the  MAGI. 

Figure  9  shows  the  deformation  configurations  for  the  aluminum  target  and  steel  projectile  at  1.0  ps 
predicted  using  the  CMAGI  and  MAGI  codes.  It  is  seen  from  the  figure  that  the  overall  responses  are 
similar  but  the  local  deformations  are  different. 

The  discrepancy  of  the  local  deformations  could  be  caused  by  the  stress  scaling  back  algorithm.  The 
MAGI  code  brings  the  stress  deviators  5„p  back  to  the  yield  surface  from  the  overshooting,  elastic  stress 
state  based  on  the  assumption, 

where  are  the  overshooting  deviatoric  stresses  evaluated  using  the  linear  elasticity;  and  are  the 
von  Mises  stress  and  the  material  yield  stress,  respectively.  It  should  be  noted  that  Eq.  (57)  is  suitable 
only  for  the  cases  that  all  deviatoric  stresses  are  increasing  in  ratio.  In  the  impact  events,  the  loading 
could  be  history-sensitive;  therefore,  the  accuracy  of  using  Eq.  (57)  is  questionable.  In  the  CMAGI,  the 
elastic-plastic  stress  calculation  is  carried  out  using  the  scaling  factor  approach  as  previously  mentioned. 
Obviously  the  present  approach  would  predict  more  accurate  results. 

The  other  possible  reason  for  the  different  local  deformations  could  be  the  stress  and  strain 
transformations  between  the  undeformed  and  deformed  configurations.  No  transformation  is  included 
in  the  MAGI  code. 
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(a)  isotropic  material  approach 


(b)  composite  laminate  approach 


Fig.  9  Aluminum  target  and  steel  projectile  deformation  at  1.0  |is 
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4.2  SPH  SOLUTION  ACCURACY  AS  A  FTJNCTION  OF  NUMBER  OF  PARTICLES 


Two  different  AS4/3501-6  composite  laminates,  [Og]  and  [0/90/45/-45]s,  were  considered  for  studying 
the  resolution  dependency  of  the  SPH  solution  as  a  function  of  number  of  particles.  Normal  impact  at 
the  center  of  the  plates  and  plane  strain  condition  were  assumed.  The  width  and  thickness  of  the  cross- 
section  of  both  composite  targets  are  4.0  cm  and  0.0965  cm,  respectively.  The  size  of  the  square  steel 
projectile  is  0.386  cm.  The  striking  speed  of  the  projectile  was  set  at  2.0  km/s.  On  account  of  the 
symmetry,  only  half  of  the  domain  was  analyzed.  For  each  impact  case,  two  different  SPH  models  were 
set  up.  One  model  consisted  of  six  particles  per  layer,  and  the  other  model  had  two  particles  per  layer. 
Of  course,  the  simulation  with  more  particles  required  more  computer  memory  and  time  than  that  with 
less  particles. 

Figure  10  presents  the  impact  results  at  3  .0  ps  for  the  [Og]  laminate  modeled  with  six  particles  per  layer. 
The  pink-colored  particles  in  Fig.  10b  are  damaged.  Delamination  is  clearly  present  between  the  2nd 
and  3rd  layers  (from  the  bottom).  The  damaged  zone  is  about  three  times  the  projectile  size.  The  results 
predicted  with  two  particles  per  layer  are  shown  in  Fig.  11.  Although  delamination  is  not  clearly  shown 
as  in  Fig.  10,  the  patterns  and  sizes  of  the  damage  predicted  by  the  two  models  are  close  . 

The  impact  results  at  3.0  ps  for  the  [0/90/45/-45]s  laminate  predicted  with  6  particles  and  2  particles  per 
layer  are  shown  in  Figs.  12  and  13,  respectively.  The  sizes  of  the  perforation  holes  are  close  for  the  two 
predictions;  however,  the  six  particles  model  results  in  more  damage  than  the  two  particles  model. 

The  above  results  imply  that  in  order  to  accurately  predict  damage,  a  composite  structure  with  different 
fiber  orientations  or  material  systems  requires  finer  particles  than  a  single  material  structure. 

4.3  3-D  STMIJLATTONS  OF  COMPOSITE  LAMINATE  PENETRATION 

All  the  WL  AS4/3501-6  composite  specimens  measured  20.32x20.32  cm  and  were  clamped  around  the 
periphery.  The  unsupported  surface  was  17.78x17.78  cm.  Each  laminar  is  0.01206  cm  thick.  Two 
different  fragment  projectiles  were  of  interest.  One  was  prismatic  with  the  dimensions  of 
0.584x0.584x0.724  cm,  and  the  other  was  cubic  with  a  side  length  of  0.970  cm. 

As  shown  in  Section  4.2,  a  SPH  model  with  two  particles  per  layer  is  not  able  to  predict  accurate 
damage  results  for  the  [0/90/45/-45]s  laminate.  For  the  WL  8-ply  laminates,  the  ratio  between  the  length 
of  the  free  surface  and  thickness  is  about  184.  Thus,  the  number  of  particles  required  to  modeled  a 
quadrant  of  the  plates  will  be  over  8.6  millions  even  though  only  two  particles  are  modeled  in  one  layer. 
This  makes  it  formiable  difficult  to  simulate  the  WL  penetration  tests.  Therefore,  it  was  our  intent  to 
conduct  the  3-D  composite  laminate  penetration  calculations  for  demonstrating  the  SPH  capability  rather 
than  for  simulating  the  real  WL  tests. 

In  the  following  3-D  calculations,  normal  impact  simulations  of  a  steel  fragment  striking  at  the  center 
of  the  AS4/3  506-1  composite  plates  were  carried  out.  The  laminate  considered  here  is  an  8-layer 
composite  with  the  in-plane  dimensions  of  4. 0x4.0  cm.  The  cubic  projectile  is  a  steel  fragment.  Each 
size  is  0.386  cm  long  (four  times  the  laminate  thickness).  No  constraints  are  imposed  on  the  composite 
target.  For  convenience,  only  the  first  quadrant  of  the  plate  was  modeled  with  one  particle  per  layer.  The 
numbers  of  particles  composing  the  projectile  and  the  composite  are  8,192  and  217,800,  respectively. 
It  should  be  noted  that  the  quasi-isotropic  laminate  considered  here  exhibits  the  bending  twist  coupling. 
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Fig.  10  8-layer  unidirectional  composite  impacted  by  a  steel  rod  at  3.0  p-S,  V=2.0  km/s 
(6  particles  per  layer) 
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(2  particles  per  layer) 


29 


Fig.  12  The  [0/90/45/-45]s  laminate  impacted  by  a  steel  rod  with  V=2.0  km/s  at  3.0  |is 
(6  particles  per  layer) 
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Fig.  13  The  [0/90/45/-45]s  laminate  impacted  by  a  steel  rod  with  V=2.0  km/s  at  3.0  ps 
(2  particles  per  layer) 


Strictly  speaking,  an  entire  plate  or  a  half  of  the  plate  with  a  special  treatment  of  the  cut  boundary  should 
be  analyzed. 

Case  1:  [0/90/45/-45]s  laminate  impacted  at  2.0  km/s 

Figure  14  presents  the  predicted  deformation  states  at  different  penetration  stages  for  the  [0/90/45/-45], 
laminate  impacted  at  2.0  km/s.  The  length  of  each  size  of  the  plate  shown  in  the  figure  is  0.6  cm.  As 
seen  in  Figs.  14c  and  d,  the  perforation  process  completed  between  1.0  ps  and  2.0  ps.  Ejectas  are  found 
in  Fig.  14d.  They  are  originated  from  the  top  layer  where  the  fibers  are  perpendicular  to  the  size  of  the 
projectile.  Figure  15  shown  the  different  views  of  the  deformation  configurations  at  3.0  ps.  The  damage 
at  this  time  is  considered  as  the  final  damage  state  although  there  was  a  small  amount  of  additional 
deformation  still  to  occur.  It  is  visually  indicated  that  the  composite  laminate  remains  in  good  shape 
except  for  the  vicinity  of  the  impact  region.  The  damage  pattern  at  the  cross-sections  of  x— 0,  y=0,  and 
x=y  are  plotted  in  Fig.  16.  The  length  of  the  cross-section  shown  in  this  figure  is  1.0  cm.  The  particles 
with  yellow  color  are  damaged  with  at  least  one  failure  mode.  As  shown  in  Fig.  16,  the  perforation  hole 
is  about  30%  larger  than  the  size  of  the  projectile.  Although  the  present  model  of  one  particle  per  layer 
can  not  exactly  predict  delamination,  the  damage  pattern  at  y=0  (i.e.  x-axis)  is  similar  to  that  reported 
by  Czamecki  [36].  In  Ref  36,  the  quasi-isotropic  32-ply  laminate  was  made  of  the  same  [0/90/45/-45] 
family.  The  damage  sustained  near  the  laminate  rear  face  was  extreme  delamination  several  times  the 
diameter  of  the  spherical  projectile. 

Case  2:  [45/-45/0/90]s  laminate  impacted  at  2.0  km/s 

The  deformation  at  3.0  ps  for  the  [45/-45/0/90],  laminate  impacted  at  2.0  km/s  is  shown  in  Fig.  17.  The 
length  of  each  size  of  the  plate  shown  in  the  figure  is  0.6  cm.  No  ejecta  is  found  in  this  case.  The  size 
of  the  perforation  is  close  to  that  in  the  [0/90/45/-45]s  laminate  impacted  by  the  same  speed.  Figure  18 
indicates  the  damage  patterns  at  the  cross-sections  of  x=0,  y=0,  and  x=y.  The  length  of  the  cross-section 
shown  in  Fig.  18  is  1.0  cm.  The  detailed  damage  in  the  vicinity  of  the  perforation  is  different  from  that 
seen  in  the  [0/90/45/-45]s  laminate.  At  x=0,  the  dominating  damage  in  the  second  and  seventh  plies  is 
much  severer  in  the  [45/-45/0/90J5  laminate  than  that  in  the  [0/90/45/-45]5  laminate.  At  y=0,  the  damage 
occurs  evenly  in  the  two  -45  degree  plies  and  the  0-degree  layer  in  this  case  instead  of  the  triangular 
damage  zone  found  in  the  three  bottom  layers  for  the  [0/90/45/-45]3  laminate  (Fig.  16b).  The  difference 
is  less  pronounced  for  the  cross-section  of  x=y. 

Case  3:  [45/-45/0/90]s  laminate  impacted  at  0.6  km/s 

Figure  19  shows  the  deformation  at  10.0  ps  for  the  [45/-45/0/90],  laminate  impacted  at  0.6  km/s.  The 
length  of  each  size  of  the  plate  shown  in  the  figure  is  1.0  cm.  It  is  interesting  to  note  the  results  in  Fig. 
19b.  A  strip  of  damage  extends  along  the  45  degree  direction  from  the  comer  of  the  perforation.  This 
damage  pattern  is  similar  to  the  C-SCAN  inspection  (Fig.  20)  made  by  the  Wright  Laboratory  [37]. 
Figure  21  presents  the  detailed  damage  at  the  cross-sections  of  x=0,  y=0,  and  x=y.  The  length  of  the 
cross-section  shown  in  the  figure  is  1.5  cm.  The  front  size  of  the  perforation  is  close  to  the  size  of  the 
projectile,  while  the  rear  size  is  larger.  In  general,  the  perforation  hole  is  smaller  than  that  caused  by  the 
striking  speed  of  2.0  km/s.  The  damage  along  the  45-degree  direction  is  clearly  indicated  in  Fig.  21c. 
It  occurs  in  the  two  bottom  layers  and  the  top  layer  and  is  several  times  the  size  of  the  projectile. 

The  penetraion  simulation  for  the  second  quadrant  of  the  plate  was  also  carried  out.  No  damage  was 
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(c)  t  =  1 .0  |j.s 


(d)  t  =  2.0|is 


Fig.  14  Penetration  at  different  stages  for  the  [0/90/45/-45]s  laminate  impacted  by 
a  steel  cube  with  V=2.0  km/s 
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Fig.  15  Deformation  at  3.0  |i.s  for  the  [()/90/45/-45]s  laminate  impacted  by  a  steel  cube  with  V=2.0km/s 


Deformation  at  3.0  |j.s  for  the  [45/-45/0/90]s  laminate  impacted  by  a  steel  cube  with  V=2.0  km/s 


=0.6  km/s 


(a)  x=0 


(b)  y=0 


(c)  x=y 


Fig.  2 1  Damage  pattern  at  the  cross-sections  of  x=0,  y=0,  and  x=:y  in  the  [45/-45/0/90]s 
laminate  at  10.0  ps  after  impact  by  a  steel  cube  at  V=  0.6  km/s 
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found  on  the  top  and  bottom  faces  except  in  the  vicinity  of  the  perforation.  This  is  consistent  with  the 
experimental  result  shown  in  Fig.  20. 

Case  4:  [45/-45/0/90]s  laminate  impacted  at  0.2  km/s 

The  deformation  at  10.0  ps  for  the  [45/-45/0/90],  laminate  impacted  at  0.2  km/s  is  presented  in  Fig.  22. 
The  length  of  each  size  of  the  plate  shown  in  the  figure  is  1.0  cm.  The  penetration  process  is  not 
completed  yet  at  this  time.  Note  again  the  strip  of  damage  that  extends  along  the  45  degree  direction 
fi-om  the  comer  of  the  perforation  (Fig.  22b).  The  width  of  the  strip  is  greater  than  that  found  in  the  case 
of  the  0.6  km/s  impact  speed.  The  final  damage  in  the  entire  quadrant  of  the  composite  plate  (at  30.0  ps) 
is  shown  in  Fig.  23.  The  perforation  hole  is  larger  in  the  x-direction  than  in  the  y-direction.  The  damage 
is  found  at  many  locations  in  both  Figs.  23a  and  b.  It  is  not  clear  if  the  prediction  is  reasonable  or  not 
for  this  case. 

Figure  24  shows  the  deformed  projectiles  after  impact  at  the  striking  speeds  of  0.2,  0.6  and  2.0  km/s. 
It  is  clearly  indicated  in  the  figure  that  the  projectile  loses  its  mass  at  2.0  km/s.  The  higher  the  impact 
speed,  the  severer  the  deformation  of  the  projectile.  The  predicted  residual  velocities  of  the  projectiles 
are  about  0.185  km/s,  0.558  km/s  and  1.885  km/s,  respectively.  The  percentages  of  the  residual 
velocities  to  the  initial  striking  speeds  are  92.5%,  93.0%  and  94.2%,  which  are  close  to  the  measurement 
by  the  Wright  Laboratory  [37]. 
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Fig.  22  Deformation  at  10.0  |j.s  for  the  [45/-45/0/90]s  laminate  impacted  by  a  steel  cube  at  V=0.2  km/s 
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Fig .  23  Damage  in  the  [45/-45/0/90]s  laminate  at  30.0  |a.s  after  impact  by  a  steel  cube  at  V=0.2  km/s 
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Fig.  24  Deformation  of  the  projectiles  after  impact  at  different  striking  speeds 


5.  CONCLUSIONS  AND  RECOMMENDATION 


A  physics-based  "first  principle"  model  for  composite  laminate  penetration  was  developed  and 
incorporated  into  the  smoothed  particle  hydrodynamics  code.  The  model  includes  various  modules 
appropriate  for  fiber  composites.  They  are  an  anisotropic  elasto-plasticity  constitutive  model,  an 
equation  of  state,  a  failure  criterion,  as  well  as  stress  and  strain  transformations  between  the  material  and 
geometric  coordinate  systems  under  large  deformation.  Due  to  the  limitation  of  the  computer  memory 
and  time,  numerical  simulations  of  3-D  composite  laminate  penetration  were  performed  with  an  intent 
for  the  SPH  capability  demonstration  rather  than  the  real  WL  test  simulation.  Our  results  clearly 
demonstrate  success  in  both  accounts:  (a)  numerical  simulation  of  composite  laminate  damage  zone 
correlated  well  with  the  known  experimental  data,  and  (b)  newly  improved  SPH  proved  to  be  a  robust 
and  viable  analytical  tool  for  predicting  response  of  non-homogeneous,  anisotropic  composite  materials 
to  intensive  dynamic  loadings. 

In  conclusion,  we  have  achieved  the  established  goal  of  this  task  where  a  fundamental  and  economical 
analytical  tool  capable  of  predicting  damage  response  of  fiber-reinforced  composites  to  impact  loadings 
is  developed.  We  submit  however,  there  is  a  great  deal  more  development  needed  before  we  can  assess 
damage  across  a  broad  spectrum  of  composite  structures.  Following  is  a  list  of  most  pertinent  topics; 

(1)  generalization  of  the  current  developed  anisotropic  plasticity  model  and  failure  criterion  for 
fiber-reinforced  composites  to  encompass  high  strain-rate  regime, 

(2)  extension  of  the  proposed  equation  of  state  for  composites  to  include  temperature  and  phase 
changes, 

(3)  development  of  a  tri-axial  ellipsoidal  smoothing  function  for  prismatic  particles,  and 

(4)  extension  of  the  current  SPH  capability  for  interface  and  boundary  condition  applications. 

The  first  two  topics  are  straight  forward.  They  heavily  rely  on  experiment.  The  last  two  needs  are 
explained  as  follows. 

To  overcome  the  computer  memory  problem  and  to  improve  the  computational  efforts,  parallel 
computation  could  be  one  of  the  possible  approaches.  However,  this  method  could  be  less  attractive 
because  the  number  of  particles  required  for  an  accurate  solution  could  be  enormous.  An  alternative 
approach  is  the  mixed  finite  element  and  SPH  method.  The  critical  regions,  such  as  the  vicinity  of  the 
projectile-target  interaction  and  the  possible  damage  zones,  are  modeled  using  the  SPH,  while  the 
remote,  non-critical  area  is  analyzed  using  the  finite  element  method.  Again,  this  approach  could  be  less 
intriguing.  The  reason  is  that  damage  zone  in  composite  laminates  could  be  large  and,  in  turn,  requires 
a  huge  number  of  particles  to  model.  The  other  reason  is  that  the  boundary  of  the  damage  zone  is 
unknown  prior  to  the  solution.  Of  consequence,  the  partition  of  the  entire  domain  for  the  finite  element 
method  and  SPH  depends  on  one's  luck.  A  more  feasible  approach,  we  believe,  is  to  develop  a  tri-axial 
ellipsoidal  smoothing  function  for  prismatic  particles  [38-39]  .  A  particle  with  a  large  aspect  ratio  of  the 
in-plane  dimension  to  thickness  can  tremendously  reduce  the  total  number  of  particles  needed  to  model 
the  composite  laminate  and  the  projectile.  A  technical  problem  of  this  approach  is  that  an  ellipsoidal 
kernel  does  not  conserve  angular  momentum.  Once  the  problem  is  solved,  the  benefits  from  this 
approach  on  computer  memory  and  time  as  well  as  labor  efforts  are  evident. 

Because  the  damage  in  a  composite  laminate  could  significantly  be  influenced  by  the  stress  waves 
reflected  from  different  types  of  boundary,  the  boundary  conditions  should  properly  be  treated  in  high 
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velocity  impacts.  In  the  current  CMAGI  code,  two  boundary  conditions  are  allowed:  free  and  symmetric. 
The  latter  is  imposed  by  introducing  ghost  particles.  To  treat  the  clamped  boundary  conditions  as 
required  by  the  WL  tests,  some  additional  technical  work  must  and  can  be  done. 
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